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Abstract 

We provide a brief overview on the application of the exterior calculus of differential forms to 

S-H . 

^H' the ab initio formulation of field theories based upon random simplicial lattices. In this framework, 
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< 



discrete analogues of the exterior derivative and the Hodge star operator are employed for the 
factorization of discrete field equations into a purely combinatorial (metric- free) part and a metric- 



dependent part. The Hodge star duality (isomorphism) is invoked to motivate the use of primal 

and dual lattices (a dual cell complex) . The natural role of Whitney forms in the construction of 
-(— > \ 

discrete Hodge star operators is stressed. 
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I. INTRODUCTION 

The need to formulate field theories on a lattice (mesh, grid) arises from two main reasons, 
which may occur simultaneously or not. First, the lattice provides a natural 'regularization' 
of divergences in lieu of renormalization techniques (l|. Such regularization does not need 
to be viewed as an ad hoc step, but instead as a natural consequence of assuming the field 
theory to be, at some fundamental level, an effective ('low'-energy) description J2j]. Second, 
the lattice provides a direct route to compute, in a non-perturbative fashion, quantities of 
interest by numerical simulations. Nontrivial domains and complex boundary conditions can 
then be easily treated as well [3[] , [4J , [5j , [6| . For these, the use of irregular ('random') lattices 
are often of interest to gain geometrical flexibility. Irregular lattices are also of interest as 
a means to provide a potentially faster convergence to the continuum limit, near- isotropic 
lattice dispersion properties, and better 'conservation' of some (e.g., long-range translational 
and rotational) symmetries |7|,[8|. In some cases, irregular lattices are useful for universality 
tests as well [9|, 10]. 



Lattice theories are typically developed by taking the counterpart continuum theory as 
starting point and then applying discretization techniques whereby derivatives are approx- 
imated by finite-differences or some constraints are enforced on the functional space of ad- 
missible solutions to be spanned by a finite set of 'basis' functions (e.g., 'Galerkin methods' 
such as spectral elements and finite elements). These discretization strategies have proved 
very useful in many settings; however, they often produce difficulties in the case of irregular 
('random') lattices. Among such difficulties are (i) numerical instabilities in marching-on- 
time algorithms (regardless of the time integration method used), (ii) convergence problems 
in algorithms relying on iterative linear solvers, and (Hi) spurious ('ghost') modes and/or 
extraneous degrees of freedom. These problems often (but not always) appear associated 
with highly skewed or obtuse lattice elements, or at the boundary between heterogeneous 
(hybrid) lattices subcomponent, comprising overlapped domains or "mesh- stitching" inter- 
faces, for example. Clearly, such difficulties put a constraint on the geometric flexibility 
that irregular lattices are intended for, and may require stringent (and computationally de- 
manding) mesh quality controls. These difficulties also impact the ability to utilize 'mesh 
refinement' strategies based on a priori error estimates. The reasons behind these difficulties 
can be traced to an inconsistent rendering of the differential calculus and degrees of freedom 



on the lattice. A rough classification of those inconsistencies is provided in Appendix D. 

The objective of this work is to provide a brief overview on the ap- 
plication of exterior calculus of differential forms to the ab initio for- 



mulation of 
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d theories on irregular simplicial ( or 'random') la; 
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In the exterior calculus framework, the lattice is treated as a cell complex (in the parlance 

n 

of algebraic topology [30[) instead of simply a collection of discrete points, and dynamic 
fields are represented by means of discrete differential forms (cochains) of various de- 



grees 



, 31], 32j. This prescription provides a basis for developing a consistent 'discrete 
calculus' on irregular lattices, and discrete analogues to partial differential equations that 
better adheres to the underlying physics. 

This topic intersects many disparate application areas. For concreteness, we use classical 
electrodynamics in 3+1 dimensions as a basic model. Although some familiarity with the 
exterior calculus of differential forms is assumed [18|], [l9(] ,[33j], [34] ,[35(], [36(], [37|],[38|], J39|] , the 



discussion is mostly kept at a tutorial level. Finally, we stress that this is a review paper 
and no claim of originality is intended. 



II. PRE-METRIC LATTICE EQUATIONS 

Let us denote the space of differential p-forms on a smooth connected manifold Q as 
A P (Q). From a geometric perspective, a differential p-form a p G A P (Q) can be viewed as 
an oriented p- dimensional density, or an object naturally associated with p-dimensional 
domains of integration U p such that the lattice contraction ('pairing') below: 



(U p ,a p )= / a 



(1) 



gives a real number (in our context) for each choice of U p [22|. On a lattice /C, U p is restricted 
to be a union of elements from the finite set of p-dimensional N p oriented lattice elements, 
which we denote r p (/C) = {a. 



p,« ) 



1, . . . , N p }. These are collective called 'p-chains'. In 
four-dimensions for example, they correspond to the possible unions of elements from the set 
of vertices (nodes) a , edges ('links') Oi, facets ('plaquettes') cr 2 , volume cells ('voxels') a 3 , 
and hypervolume cells 04, for p = 1, . . . , 4, respectively. In the discrete setting, the degrees 
of freedom are reduced to the set of pairings (CD) on each one of the lattice elements. 



On the lattice, the pairing above can be understood as a map W : A P (Q) — > T p (/C) such 
that 

W{a p ) = ((T Pti ,o^)= f a p (2) 



defines its action on the basis of p-chains. Note that we use T p (/C) to denote the space 
dual to r p (/C), i.e. the space p-cochains. The latter can be viewed as the space of 'discrete 
differential forms'. Because of this, and with some abuse of language, we use the terminology 
'differential forms' and 'cochains' interchangeably to denote the same objects in what follows. 



The map TZ P is called the de Rham map 22 1. 



The basic differential operator of exterior calculus is the exterior derivative d, applicable 
to any number of dimensions. The discretization of d on a general irregular lattice can be 



effected by a straightforward application of the generalized Stokes' theorem 22] 



/ da p = J a p (3) 

with p = 0, . . . , 3 in n = 4. In the above, d is the boundary operator, which simply maps a 
p-dimensional lattice element to the set of (p— l)-dimensional lattice elements that comprise 
its boundary, preserving orientation. This theorem sets d as the formal adjoint of d in terms 
of the pairing given in (jTJ, that is (cr p+ i, da p ) = (da p+ i, a p ). Computationally, the boundary 



operator can be implemented by means of incidence matrices 22J, 29), 40( such that 



3 

where the indices i and j run over all (p+1)- and p-dimensional lattice elements, respectively. 
The incidence matrix entries are such that Cf, G { — 1,0, 1} for all p, with sign determined 
by the relative orientation of lattice elements i and j. The restriction to this set of integer 
values reflects the 'metric-free' nature of the exterior derivative: only information about 
element connectivity, that is, the combinatorial aspects of the lattice, is involved here. It 
turns out that the metric is fully encoded by Hodge star operators, the discretization of 
which will be discussed further down below. 
Using eqs. (J3J) and (@|, one can write 



da p 



E^7 aP ( 5 ) 

J Or, A 



for all i, so that the derivative operation is replaced by a proper sum over j. On the lattice, 
the nilpotency of the operators dod = dod = [41] is recovered by the constraint [22 1 



E c s hlc & = ° 



(6) 



for all i and j. 



ILL EXAMPLE: LATTICE ELECTRODYNAMICS 



We write Maxwell's equations in a four- dimensional Lorentzian manifold Q as 34] 



dF = 



dG = *J 



(7) 
(8) 



where d is the four- dimensional exterior derivative, F and G are the so-called Faraday and 
Maxwell 2-forms, respectively, and *J is the charge- current density 3-form. The Hodge star 
operator * is an isomorphism that maps p-forms to (4 — p)-forms, and more generally p 



forms to (n — p) forms in a n-dimensiona 
the metric of ft [231 > M > W > M > M > Wl > 



manifold, and, as mentioned before, depends on 



44j, 45|. The above equations are complemented 



by the relation G = *F, which indicates that F and G are 'Hodge duals' of each other. 



Primal and dual lattices 



Since F and G are 2-forms, they should be discretized as 2-cochains residing on plaquettes 
(2-chains) of the 4-dimensional lattice; however, it is important to recognize that these two 
forms are of different types: F is a 'ordinary' (or 'non-twisted') differential form, whereas G 
(as well as *J") is a 'twisted' (or 'odd') differential form [46]. The basic difference here has 
to do with orientation: ordinary forms have internal orientation whereas twisted forms have 



external orientation 20], 22] , [46] , [47j , [48] . These two types of orientations exhibit different 
symmetries under reflection, a distinction akin to that between proper (or polar) tensors and 
pseudo (or axial) tensors. Only twisted forms admit integration in non-orientable manifolds. 
These two types of forms are associated with two distinct 'cell complexes' (lattices), each 
one inheriting the corresponding orientation: the ordinary form F is associated with the 
set of plaquettes T 2 on the 'ordinary cell complex' /C, thus belonging to T 2 (/C), while the 



twisted forms G and *J are associated with the set of plaquettes T 2 on the 'twisted cell 



complex' /C [22J , [27| , |48j , 49], thus belonging to T 2 (/C). Consequently, we also have two sets 



of incidence matrices Cf, and Cf,, one for each lattice. It is convenient to denote /C as the 
'primal lattice' and K, as the 'dual lattice' [221 ] . 

As detailed further below, these two lattices become intertwined by the Hodge duality F = 
*G. The need for dual lattices can also be motivated from a purely combinatorial standpoint 



(as a means to preserve key topological properties from the continuum theory) 



241 or from 



a strictly computational 



for example) [501]. [5 1|. [52]. 



standpoint (to provide higher-order convergence to the continuum, 



B. 3+1 theory 

At this point, it is suitable to degeometrize time and treat it simply as a parameter. 
This corresponds to the majority of low-energy applications involving Maxwell's equations, 
in which one is interested in predicting the field evolution along different spatial slices for 
a given set of initial and boundary conditions. In this case, we still use the symbols K, and 
K, for the primal and dual lattices, but they now refer to three-dimensional spatial lattices. 
Similarly, Q now refers to a three-dimensional Euclidean manifold . In such a 3+1 setting, 
one can decompose F and G as 

F = EAdt + B (9) 

G = D-H Adt (10) 

and the source density as 

*J = -JAdt + p (11) 

where A is the wedge product, E and H are the electric intensity and magnetic intensity 
1-forms on 1+ and r\ respectively, D and B are the electric flux and magnetic flux 2-forms 
on T2 and 1^ respectively, J is the electric current density 2-form on 1^ , and p is the 
electric charge density 3-form on T 3 (corresponding assignments for the 2+1 and 1+1 cases 
are provided in |32J). As a result, Maxwell's equations reduce to 

dE = -d t B (12) 

dH = d t D + J (13) 



representing Faraday's and Ampere's law, respectively. Here, d stands for the 3- dimensional 
spatial exterior derivative. Note that both eqs. f JT2|) and fflBl are metric-free. They are 
supplemented by Hodge star relations given by 



D = * f E 



H = -k^-iB 



(14) 
(15) 



now involving two Hodge star maps in three-dimensional space: * e : A x (f2) — > A 2 (Q) and 
* M -i : A 2 (Q) — > A l (Q). On the lattice, we have the corresponding discrete counterparts: 
[* e ] : T 1 (/C) — ¥ T 2 (/C) and [* M -i] : T 2 (/C) — > T 1 (/C). The subscripts e and fi in * e and * M -i serve 
to indicate that these operators also incorporate macroscopic constitutive material properties 
through the local permittivity and permeability values 53| (we assume dispersionless media 
for simplicity). In Riemannian manifolds (and in particular, Euclidean space ) an d reciprocal 
media, these two Hodge star operators are symmetric and positive-definite [54 1. 

In what follows, we employ the following short-hand notation for cochains: (<Ti t i, E) = Ei, 
{o\,i,H) = Hi, (a 2 ,i,D) = Di, (a 2 ,i,B) = B h (a 2 ,i,J) = Ji, and (a 3ji ,p) = p h where the 
indices run over the respective basis of p-chains in either K, or K., p — 1, 2, 3. With the 
exception of Appendix A, we restrict ourselves to the 3+1 setting throughout the remainder 
of this paper. 



IV. CASTING THE METRIC ON A LATTICE 



A. Whitney forms 



The Whitney map W : T p (/C) — > A P (Q) is the right-inverse of the de 
Rham map (T5]), that is, TZ o W = X, where X is the identity operator. In 
simplicial lattices, this moronism can be constructed using the so-called Whitney 



forms (l5j,|22|,|36|,|43|,|55|,|56 



J57l| . J58| . J59l| . |60| . J6l| which are basic interpolants from 
cochains to differential forms 33] (other interpolants are also possible 62|,63[|). By def- 
inition, all cell elements of a simplicial lattice are simplices, i.e., cells whose boundaries are 
the union of a minimal number of lower-dimensional cells. In other words, O-simplices are 
nodes, 1-simplices are links, 2-simplices are triangles, 3-simplices are tetrahedra, and so on. 



Note that if the primal lattice is simplicial, the dual lattice is not 31] . For a p-simplex <7 Pi ;, 



the (lowest-order) Whitney form is given by 



CJ 



P l a P,i] = P ! /X—iy^jdKo A d\,i ' ' ' d\i,j-i A dXij+x ■■■d\ 

j=0 



i,P 



(16) 



where Xi t j, j = 0, . . . ,p, are the barycentric coordinates associated to a v ^. In the case of a 
0-simplex (node), f|T6|) reduces to w°[cro,i] = Aj. 

From its definition, it is clear that Whitney forms have compact support. Among its 
important structural properties are: 



(a p ,i,uj p [a p j}) 






'I.I 



(17) 



where 5ij is the Kronecker delta, which is simply a restatement of TZ o W = X, and 



(18) 



where d T is the coboundary operator [56], consistent with the generalized Stokes' theorem. 



Further structural properties are provided in [57j,[58|]. Higher-order version of Whitney 
forms also exist 59( , 60| ■ The key result W o 1Z — y X holds in the limit of zero lattice spac- 
ing. This is discussed, together with other related convergence results in various contexts, 

in luty ,y ,y,y ,|67J,|68|. 

Using the short-hand a; p [cr P) j] = uf, we can write the following expansions for E and B 
in a irregular simplicial lattice, in terms of its cochain representations: 



i 



(19) 



(20) 



where the sums run over all primal lattice edges and faces, respectively. 

One could argue that Whitney forms are continuum objects that should have no funda- 
mental place on a truly discrete theory. In our view, this is only partially true. In many 
applications (see, for example, the discussion on space-charge effects below), it is less natural 
to consider the lattice as endowed with some a priori discrete metric structure than it is to 
consider it instead as embedded in an underlying continuum (say, Euclidean) manifold with 
metric and hence inheriting all metric properties from it. In the latter case, Whitney forms 
provide the standard route to incorporate metric information into the discrete Hodge star 
operators, as described next. 
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B. Discrete Hodge star operator 

In a source-free media, we can write the Hamiltonian as 

H = \ f (EAD + H AB)= I {E A * t E + *„_i.B A B) (21) 

Using eqs. f|T9|) and (1201) . the lattice Hamiltonian assumes the expected quadratic form: 

n = E E ^ [*J« *y + E E ^ tv- 1 !* ^' ( 22 ) 

» i i 3 

where we immediately identify the symmetric positive definite matrices 

[*J«= [ul^Wj (23) 

[*^]<i= /"(•a.-iw^AwJ (24) 



as the discrete realization of the Hodge star operator (s) on a simplicial lattice 23|], 69( so 

that 

A = EN;;^- ( 25 ) 

^ = Eh- 1 ^;- ( 26 ) 

From the above, the Hamiltonian can be also expressed as 

u = J2 EiDi+ J2 HiBi ( 27 ) 

i i 

C. Symplectic structure and dynamic degrees of freedom 

The Hodge star matrices [* e ] and [vr^-i] have different sizes. The number of elements in 
[* e ] is equal to Ni x i\q, whereas the number of elements in f*^ 1 ] is equal to N 2 x N 2 . In 
other words, O(-E) = 6(-D) ^ Q(B) = Q(H), where denotes the number of (discrete) 
degrees of freedom in the corresponding field. 

One important property of a Hamiltonian system is its symplectic character, associated 
with area preservation in phase space. The symplectic character of the Hamiltonian in 
principle would require a canonical pair such as E, B to have identical number of degrees of 
freedom. This apparent contradiction can be explained by the fact that Maxwell's equations 
ffi~2]) and ffT3|) can be thought as a constrained dynamic system (by the divergence conditions) 



so that, even though Q(E) ^ O(-B), we still have Q d (E) = Q d (B) } where <d d denotes the 
number of dynamic degrees of freedom. This is discussed further below in Section VI, in 
connection with the discrete Hodge decomposition on a lattice. 

V. SEMI-DISCRETE EQUATIONS 

A. Local and ultra-local lattice coupling 

By using a contraction in the form of (|2J) on both sides of (1T21) with every face o"2j of /C, 
and using the fact that (cr 2 j, of) = {cf\,j ■,(*}}) = Sij from (fTTj) . we get 

(a 2j , d t B) = d t J2 B i (*2j, w?) = d tBj (28) 

i 

and 

K,, dE) = (da 2>j , E) = Y J E,Y, tfk (°i,k, "]) = J2 C h E > ( 29 ) 

i k i 

so that 

-d t B l = J2cl 3 E 3 (30) 

3 

where the index i runs over all faces of the primal lattice. On the dual lattice /C, we can 
similarly contract both sides of eq. ( fT3j) with every dual face c^j to get 

$A = 5^CjJ7i (31) 

i 

where now the index i runs over all faces of the dual lattice. Using eqs. (125]) and (126]) and 

the fact that, in three-dimensions C}, = Cj { [22J (up to possible boundary terms ignored 

here), we can write the last equation in terms of primal lattice quantities as 

dt £>J« Ej = J2 C h YfrrAikBk (32) 

3 J k 

or, by using the inverse Hodge star matrix [* e ]^ , as 

d t E i = ^2T ij B j (33) 

j 

with 

T^EE^^M'i (34) 



it z 
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The matrix [T] can be viewed as the discrete realization, for p = 2, of the codifferential 



operator 5 = (— l) p * 1 d* that maps p-forms to (n — p)-forms 351 ] . 



n 

Since the continuum operators * e and ^-i are local [46] and, as seen, Whitney forms 
(TTB|) have local support, it follows that the matrices [* 6 ] and [*^-i] are sparse, indicative of 



£ 



an ultra- local coupling (in the terminology of [70]). In contrast, the numerical inverse [* e ] _1 
used in eq. (|3"4"|) is, in general, not sparse so that the field coupling between distant elements 
is nonzero. The lack of sparsity is a potential bottleneck in practical simulations. However, 
because the coupling strength in this case decays exponentially |29|,j44j, we can still say 
(using again the terminology of |70j) that the resulting discrete operator encoded by the 
matrix in (j34p is local. In practical terms, the exponential decay allows one to set a cutoff 
on the nonzero elements of [* e ], based on element magnitudes or on the sparsity pattern 
of the original matrix [* e ], to build a sparse approximate inverse for [* e ] and hence recover 
back an ultra-local representation for -k' 1 [29j.[7l|. The sparsity pattern of [* £ ] encodes the 
nearest-neighbor edge information of the mesh and, consequently, the sparsity pattern of 
{k 6 ] k likewise encodes successive '/c-level' neighbors. The latter sparsity patterns can be used 



to build, quite efficiently, sparse approximations for [* e ] _1 , as detailed in [29[. Once such 
sparse representations are obtained, eqs. f|30|) and ( 133]) can be used in tandem to construct a 
marching-on-time algorithm (see Section IX. A ahead, for example) with a sparse structure 
and hence amenable for large-scale problems. 



B. Barycentric dual and barycentric decomposition lattices 

An alternative approach, aimed at constructing a sparse discrete Hodge star for -k' 1 
directly from the dual lattice geometry is described in 27J, based on earlier ideas exposed 



m 



241 ]. 72] . This approach is based on the fact that both primal /C and dual K. lattices can 



be decomposed into a third (underlying) lattice /C by means of a barycentric decomposition, 
see 24|. The dual lattice K in this case is called the barycentric dual lattice J27J,[72j and 



the underlying lattice /C is called the barycentric decomposition lattice. Importantly, /C is 
simplicial and hence admits Whitney forms built on it using f lT6|) . Whitney forms on K, can 
be used as building blocks to construct (dual) Whitney forms on the (non-simplicial) /C, and 
from that, a sparse inverse discrete Hodge star [^T 1 ] using integrals akin to (I2"3"j) and (12"4"|) . 
An explicit derivation of such dual lattice Whitney forms is provided in 73J|. Furthermore, a 



11 



recent comprehensive survey of this and other approaches based on dual lattices to construct 
discrete sparse inverse Hodge stars is provided in |74| . 

The barycentric dual lattice has the important property below associated with Whitney 
forms: 

(a(n- P ),i,*u p [(7p,j}) = *u p [a pJ ] = 6ij (35) 



where * stands for the spatial Hodge star operator (distilled from constitutive material 

properties), and a( n - p ) t i is the dual element to <r Pji on the barycentric dual lattice. The 

operator * is such that 

[ u p A*u p = [ \co\ 2 dv (36) 

Jn Jn 

where \u\ 2 is the two- norm of u p and dv is the volume element. 

The identity (|35|) plays the role of structural property (ITT)) , on the dual lattice side. 
We stress that identity (135]) is a distinctively characteristic feature of the barycentric dual 
lattice not shared by other geometrical constructions for the dual lattice. In other words, 
compatibility with Whitney forms via ( 155|) naturally forces one to choose the dual lattice to 
be the barycentric dual. 

From the above, one can also define a (Hodge) duality operator directly on the space of 
chains, that is *k '■ T p (/C) !->■ r n _ p (/C) with *K{vp,i) = &(n-p),i an d *x : Tp(^) ^ r„_„(/C) 



with *ft-(cr Pi j) = a(n-p) t i, so that *k*k — ^k^k = 1- This construction is detailed in [24]. 



C. Galerkin duality 

Even though we have chosen to assign E and B to the primal (simplicial) lattice, and 
consequently D, H, J, and p to the dual (non-simplicial) lattice, the reverse is equally 
possible. In this case, the fields D, H become associated to a simplicial lattice and hence 
can be expressed in terms of Whitney forms; the expressions dual to f|T9|) and (1201) are now 

H = J2^ ul (37) 

i 

D = Y,D^ 2 (38) 

i 

with sums running over primal edges and primal faces, respectively, and where 

Ei = Yl^-^Dj (39) 

3 
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Bi = ][>„]<,-ffj (40) 

3 

with 

K% = [ {**-«*>!) A U] (41) 

./Q 

W<i = / u l A V"j ( 42 ) 

and the two Hodge star maps now used are such that, in the continuum, -*-" 1 : A 2 (Q) — > A l (Q) 
and *„ : A 1 ^) -* A 2 (fi), and, on the lattice, \k~ x ] : T 2 (/C) -> r 1 ^) and [* M ] : T^/C) -» 
T 2 (/C). This alternate choice entails a duality between these two formulations, dubbed 
'Galerkin duality'. This is explored in more detail in [44] . 



VI. DISCRETE HODGE DECOMPOSITION AND EULER'S FORMULA 

For any p-form a p , we can write 

a P = dC' 1 + 8(3 P+1 + x p , (43) 



where \ v 1S a harmonic form 31]. This Hodge decomposition is unique. In the particular 
case of the 1-form E, we have 

E = d<P + 8 A + X , (44) 

where is a 0-form and A is a 2-form, with dcj) representing the static field, 8 A the dynamic 
field, and \ the harmonic field component (if any). In a contractible domain, % is identically 
zero and the Hodge decomposition simplifies to 

E = d(j) + 8 A. (45) 

more usually known as Helmholtz decomposition in three-dimensions. 

In the discrete setting, the degrees of freedom of are associated to the nodes of the 
primal lattice. Likewise, the degrees of freedom of A are associated to the facets of the 
primal lattice. Consequently, we have from ( 145]) that 

Q d {E) = N^-N^ 

= {N E - N b E ) - (N v - N^) 

= N E - N v , (46) 

13 



where Ny is the number of primal nodes, Np the number of primal edges, and Np the 
number of primal facets, with superscript b standing for boundary (fixed) elements and h 
for interior (free) elements. 

On the other hand, once we identify the lattice as a network of (in general) polyhedra, 



we can apply Euler's polyhedron formula on the primal lattice to obtain 44] 



N v - N E = 1 - N F + N P , (47) 

where Np represents the number of volume cells comprising the primal lattice. A similar 
Euler's polyhedron formula applies to the (closed, two-dimensional) boundary of the primal 
lattice 

N v - N b E = 2 - N b F , (48) 

Combining Eq. (1471) and (1481) . we have 

(N E - N b E ) - (N v - N b v ) = (N F - N b F ) - (N P - 1) . (49) 

From the Hodge decomposition fj45l) . we see that Q d (E) is 

Q d (E) = N% - Ny 1 

= (N E - N b E ) - (N v - N b v ) . (50) 

Note that the divergence free condition dB = produces one constraint on the 2-form B 
for each volume element. This constraint also span the whole lattice boundary. The total 
number of the constrains for B is therefore (N p — 1) . Consequently, we have 

Q d (B) = N™ - {N P - 1) 

= (N F - N b F ) - (N P - 1) (51) 

so that 

Q d (B) = B d (E) . (52) 

This discussion can be generalized to lattices on non-contractible domains with any number 



of holes (genus), where the identity Q d (B) = Q d (E) is also satisfied [3jJ- Moreover, from 
Hodge star isomorphism, we have Q d (D) = Q d (E) and Q d (H) = Q d (B). 

In general, we can trace a direct correspondence between quantities in the Euler's poly- 
hedron formula to the quantities in the Hodge decomposition formula. For example, each 

14 



term in the two-dimensional Euler's formula Ne = Ny + (Np — 1) + g is associated to a 
corresponding term in E = dip + 5 A + x; that is, the number of edges Ne corresponds to the 
dimension of the space of lattice 1-forms E, which is the sum of the number of nodes Ny (di- 
mension of the space of discrete 0-forms </>), the number of faces (Np — 1) (dimension of the 
space of discrete 2- forms A), and the number of holes g (dimension of the space of harmonic 



forms x). A similar correspondence can be traced on a three-dimensional lattice [3JJ. This 
correspondence provides a physical picture to Euler's formula and a geometric interpretation 
to the Hodge decomposition. 
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APPENDIX A: Differential forms and lattice fermions 

Differential p-forms can be viewed as antisymmetric covariant tensor fields on rank p. 
Therefore, the ingredients discussed above are app licable to any antisymmetric tensor field 
theory, including (pure) non-Abelian theories [72J. However, for (Dirac) fermion fields the 
situation is different and, at first, it would seem unclear how differential forms could be used 
to de SCrlb e spinors. Nevertheless, a US e M coition cae indeed be e stabllS ned Q.Q.fl 
To briefly address this point, let us consider next the lattice transcription of the (one-flavor) 
Dirac equation. In this Appendix, we work on Euclidean spacetime with h = c = 1 and 
adopt the repeated index summation convention with /i, v as coordinate indices, where x is 
a point in four-dimensional space. 

It is well known that fermion fields defy a lattice description with local coupling that 
gives the correct energy spectrum in the limit of zero lattice spacing and the correct chiral 
invariance [761 ] . This is formally stated by the no-go theorem of Nielsen-Ninomiya [77J and 



is associated to the well-known 'fermion-doubling' problem 78j. A perhaps less known fact 
is that it is possible to arrive at a 'geometrical' interpretation of the source of this difficulty 
by considering the 'generalization' of the Dirac equation (7 M <9 M + m)ip(x) = given by the 
Dirac-Kahler equation 

(d - 5)V(x) = -m*(x) (53) 
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The square of the Dirac-Kahler operator can be viewed as the counterpart of the Dirac 
operator in the sense that 

(d - 8) 2 = -(dS + 8d) = -□ (54) 

recovers the Laplacian operator in the same fashion as the Dirac operator squared does, that 
is ('~f fl d f j,) 2 = — <9 M <9 M = — □, where 7 M represents Euclidean gamma matrices. 

The Dirac-Kahler equation admits a direct transcription on the lattice because both the 
exterior derivative d and the codifferential 5 can be simply replaced by its lattice analogues, 
as discussed before. However, for the Dirac equation the analogy has to further involve the 
relationship between the 4-component spinor field ip and the object \l/. This relationship 



QF 



was first established in [16|],ll7| for hypercubic lattices and later extended to non-hypercubic 
lattices in 



10(,[79|. The analysis of 16[ and (l7| has shown that \l/ can be represented by a 



16-component complex-valued inhomogeneous differential form: 

4 

tf(x) = ^V(x) (55) 

where a°(x) is a (1-component) scalar function of position or 0-form, a l (x) = a 1 Ax)dx^ is a 
(4-component) 1-form, and likewise for p = 2,3,4 representing 2-, 3-, and 4-forms with 6-, 
4-, and 1-components respectively. By employing the following Clifford algebra product 

dx" V dx u = g" v + dx» A dx u (56) 

as using the anti-commutative property of the exterior product A, we have 

dx^ V dx u + dx v V dx" = 2g> MU (57) 

which exactly matches the anticommutator result of the 7 M matrices, 7^7^ + 7^7^ = 2g fiu . 
This suggests that dx^ plays the role of the 7^ matrix in the space of inhomogeneous differ- 
ential forms with Clifford product [80], that is 

7^ i-> dx* 1 V d,, (58) 

keeping in mind that while 7 M 9 M acts on spinors, whereas dx^ V d^ = (d — 5) acts on 
inhomogeneous differential forms. This analysis leads to a 'geometrical' interpretation of 



the popular Kogut-Susskind staggered lattice fermions |8jJ,[82| because the latter can be 



made identical to lattice Dirac-Kahler fermions after a simple relabeling of variables 17]. 
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The 16-component object \1/ can be viewed as a 4 x 4 matrix that produces a four-fold de- 
generacy with respect to the Dirac equation for ip. This degeneracy is actually not a problem 
in the continuum because there is a well-defined procedure to extract the 4-components of ip 
from those of \l/ [l6j,[17| whereby the 16 scalar equations encoded by (|53p all reduce to the 
same copy of the four equations encoded by the standard Dirac equation. This procedure 



is performed by a set of 'projection operators' that form a group 16j, 83j. On the lattice, 



however, the operators d and d, as well as * (which plays a role on the space of inhomoge- 



neous differential forms \P analogous to that of 7 5 on the space of spinors ip 84j), behave 
in such a way that their action leads to lattice translations. This is because cochains with 
different p necessarily live on different lattice elements and also because * is a map between 
different lattice elements. As a consequence, the product operation of such 'group' is not 
closed anymore. This nonclosure also stems from the fact that the lattice operators d and 5 



do not satisfy Leibnitz's rule [80f| . Because of this, the degeneracy of the Dirac equation on 
the lattice is present at a more fundamental level and is harder to extricate using the Dirac- 
Kahler description than the analogous degeneracy in the continuum. In this regard, a new 
approach to identify the extraneous degrees of freedom away from the continuum was re- 



cently described in 



In addition, a split-operator approach to solve Dirac equation based 



on the methods of characteristics that purports to avoid fermion doubling while maintaining 



chiral symmetry on the lattice was very recently put forth in [86j. This approach preserves 
the linearity of the dispersion relation by a splitting of the original problem into a series of 
one- dimensional problems and the use of a upwind scheme with a Courant-Friedrichs-Lewy 
(CFL) number equal to one, which provides an exact time-evolution (i.e. with no numer- 
ical dispersion effects) along each reduced one-dimensional problem. The main (practical) 
obstacle in this case is the need to use very small lattice elements. 

Needless to say, the topic of lattice fermions is vast and we cannot do much justice to it 
here; we have only focused on the more germane aspects to main theme of this overview. 

APPENDIX B: Absorbing boundary conditions 

In many wave scattering simulations, the presence of long-range interactions with slow 
(algebraic) decay, together with practical limitations in computer memory resources, implies 
that open-space problems necessitate the use of special techniques to suppress finite volume 
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effects and emulate, for example, the Sommerfeld radiation condition at infinity. Perfectly 
matched layers (PML) are absorbing boundary conditions commonly used for this pur- 
87J,J88|, 89j,J90|. In the continuum limit, the PML provides a reflect ionless absorption 



pose 



of outgoing waves, in such a way that when the PML is used to truncate a computational 
lattice, finite volume effects such as spurious reflections from the outer boundary are ex- 
ponentially suppressed. When first introduced in the literature [87| . the PML relied upon 
the use of matched artificial electric and magnetic conductivities in Maxwell's equations 
and of a splitting of each vector field component into two subcomponents. Because of this, 
the resulting fields inside the PML layer are rendered 'non-Maxwellian'. The PML concept 
was later shown to be equivalent in the Fourier domain (d t — > —iu>) to a complex coordi- 
nate stretching of the coordinate space (or an analytic continuation to a complex-valued 



coordinate space) [8J 



, 90] and, as such, applicable to any linear wave phenomena. 
Inside the PML, the (local) spatial coordinate ( along the outward normal direction to 
each lattice boundary point is complexified as 

< 



C^C 



sdCW 



(59) 



where s^ is the so-called complex stretching variable written as s^((,u) = <\(C) + ^c(C)/ w 
with <2£ > 1 and Q^ > (profile functions). The first inequality ensures that evanescent 
waves will have a faster exponential decay in the PML region, and the second inequality 
ensures that propagating waves will decay exponentially along ( inside the PML. As opposed 
to some other lattice truncation techniques, the PML preserves the locality of the underlying 
differential operators and hence retains the sparsity of the formulation. 

For Maxwell's equations, the PML can also be effected by means of artificial material ten- 
sors (Maxwellian PML) 9l| . In three-dimensions, the Maxwellian PML can be represented 
as a media with anisotropic permittivity and permeability tensors exhibiting stratification 
along the normal to the boundary S that parametrizes the lattice truncation boundary. 
The PML tensors properties depend on the local geometry via the two principal curvatures 
of S [92[ , [93] , [94{ . The boundary surface 5* is assumed (constructed) as doubly differen- 
tiable with non-negative radii of curvature, otherwise dynamic instabilities ensue during a 
marching-on-time evolution 95| . 

From (|59|) . the PML also admits a straightforward interpretation as a complexification of 



the metric 



38(,J96|. As a result, the use of differential forms readily unifies the Maxwellian 
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and non-Maxwellian PML formulations because the metric is explicitly factored out into the 
Hodge star operators — any transformation the metric corresponds, dually, to a transforma- 
tion on the Hodge star operators that can be mimicked by modified constitutive relations [371 ] . 
In the differential forms framework, the PML is obtained by a mapping on the Hodge star 
operators: * e — > * e and *„-i — > *„-i induced by the complexification of the metric. The 
resulting differential forms inside the PML, E,D,H,B therefore obey 'modified' Hodge re- 



lations D = * e E and B = -k^-iH, but identical pre-metric equations 



flU) and (H3J). In other 



words, (fT2|) and ( TT3j) are invariant under the transformation (1591) [38J.J96]. 

APPENDIX C: Implementation of space charge effects 

In many applications related to plasma physics or electronic devices, it is necessary to 
include space charges (uncompensated charge effects) into lattice models of macroscopic 
Maxwell's equations. This is typically done by representing the charged plasma media using 
particle-in-cell (PIC) methods that track the individual particles on the lattice [97J], 98|, 99]. 
The field/charge interaction is then modeled by (i) interpolating lattice fields (cochains) to 
particle positions (gather step), (ii) advancing particle positions and velocities in time using 
equations of motion, and (in) interpolating back charge densities and currents onto the lat- 
tice as cochains (scatter step). In general, the 'particles' do not need to be actual individual 
particles, but can be a collection thereof ('macro-particles'). To put it simply, incorporation 
of space charges requires two extra steps during the field update in any marching-on-time 
algorithm, which transfer information from the instantaneous field distribution to the par- 
ticle kinematic update and vice-versa. Conventionally, this information transfer relies on 
spatial interpolations that often violate the charge continuity equation and, as a result, lead 
to spurious charge deposition on the lattice nodes. On regular lattices, this problem can be 
corrected, for example, using approaches that either subtract a static solution (charges) from 
the electric field solution (Boris/DADI correction) or directly su btra ct the residual error on 



the Gauss law (Langdon-Marder correction) at each time step [100]. On irregular lattices, 
additional degrees of freedom can be a dded as coupled elliptical constraints to produce a 
augmented Lagrange multiplier system [101]. All these approaches necessitate changes on 
the original equations, while still allowing for small violations on charge conservation. In 
contrast, Whitney forms provide a direct route to construct gather and scatter steps that 
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satisfy charge conservation exactly even on unstructured lattices [1021 ]. [1031 ] . as explained 
next. To conform to the vast majority of the plasma and electronic devices literature, we 
once more restrict ourselves here to the 3+1 setting (although a four-dimensional analysis 
in Minkowski space would have provided a more succinct discussion). 

For the gather step, Whitney forms can be used to directly compute (interpolate) the 
fields at any location from the knowledge of its cochain values, such as in (TT9|) and ( 120]) 
for example. For the scatter step, charge movement can be modeled as the Hodge-dual 
of the current 2-form J, that is, as the 1-form -kj which can be expanded in terms of 
Whitney 1-forms on the primal lattice. Here, * represents again the spatial Hodge star 
in three-dimensions distilled from macroscopic constitutive properties. The Hodge-dual 
current associated to an individual point charge can be expressed as *J = qv , where q is 
the charge value, v is the associated velocity vector, and b is the 'flat' operator or index- 
lowering canonical isomorphism that maps a vector to a 1-form, given by the Euclidean 
metric. Similarly, point charges can be encoded as the Hodge-dual of the charge density 
3-form p, that is, as the 0-form *p, which can be expanded in terms of Whitney 0-forms 
on the primal lattice. These two Whitney maps are linked in such a way that the rate of 
change on the value of the 0-cochain representing -kp at a node is associated to the presence 
of a 1-cochain representing *J along the edges that touch that particular node, leading to 
exact charge conservation at the discrete level. To show this, consider for simplicity the 
two-dimensional case of a point charge q moving from point x^ s ' to point x^> during a time 
interval r inside a triangular cell with nodes <7o,o, o"o,i> an d ^0,2, or simply 0, 1, and 2. At 
any point x inside this cell, the 0-form *p can be scattered to these three adjacent nodes via 

3 

*p = qJ2( x > UJ i) UJ i ( 6 °) 

i=i 

where we are again using the short-hand cu [o"o i j] = cu°, and the brackets represent the pairing 

expressed by ([1]). In this case, p = and the pairing integral in ([1]) reduces to a function 

evaluation at a point. Since Whitney 0-forms are equal to the barycentric coordinates 

associated of a given node, that is (x, w°) = Aj(x), we have the scattered charge q\l = 

q\i(x^) on node i for a charge q at x^ 8 \ and, similarly, the scattered charge q\{ on node i 

for a charge q at x^\ The rate of scattered charge variation on a given node i is therefore 

equal to q{\{ — Af), where q = qjr. 

During r, the particle travels through a path I from x^ to x^\ and the corresponding 
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*J can be expanded as a sum of Whitney 1- forms tot associated to the three adjacent edges 
i] = 01,12,20, that is 

*j=?£^,4}4 (6i) 

ij 

The coefficients (i,coh) represent the (oriented) current flow along the associated oriented 
edge, that is, the cochain representation of *J along edge ij. Using (TToT) . the sum of the 
total current magnitude scattered along edges 01 and 20 that flows into node is therefore 

q (- (I,u4) + (^ui» = q /(-o4 + o4) (62) 

Using tob, = XidXj — \jd\ and Ai + A 2 + A 3 = 1, the above reduces to 

q [d\ = q(\ f - A*) (63) 

Ji 

which exactly matches the rate of scattered charge variation on node obtained before. It 

is clear that similar equalities hold for nodes 1 and 2. More fundamentally, these equalities 

are a direct consequence of the structural property ( fT8l) . 



APPENDIX D: Classification of inconsistencies in nai've discretizations 

We provide below a rough classification scheme of inconsistencies arising from naive 
discretizations of the differential calculus on irregular lattices. 

(a) Pre-metric inconsistencies of first kind: We call pre-metric inconsistencies of the first 
kind those that are related to the primal or dual lattices taken as separate objects and 
that occur when the discretization violates one or more properties of the continuum theory 
that is invariant under homeomorphisms — for example, conservations laws that relate a 
quantity on a region S with an associated quantity on the boundary of the region, dS (a 
topological invariant). Perhaps the most illustrative example is violation of 'divergence- free' 
conditions caused by improper construction of incidence matrices, whereby the nilpotency 
of the (adjoint) boundary operator, d o d = 0, is not observed. This implies, in a dual 



fashion, that the identity d 2 = is violated 22] . Stated in another way, the exact sequence 



property of the underlying de Rham differential complex is violated 104J . In practical 
terms, this leads to the appearance spurious charges and/or spurious ('ghost') modes. As 
the classification suggests, these properties are not related to metric aspects of the lattice, 
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but only to its "topological aspects" that is, on how discrete calculus operators are defined 
vis-a-vis the lattice element connectivity. In more mathematical terms, one can say that the 
structure of the (co)homology groups of the continuum manifold is not correctly captured 
by the cell complex (lattice). We stress again that, given any dual lattice construction, 
pre-metric inconsistencies of the first kind are associated to the primal or dual lattice taken 
separately, and not necessarily on how they intertwine. 

(b) Pre-metric inconsistencies of second kind: The second type of pre-metric inconsistency 
is associated to the breaking of some discrete symmetry of the Lagrangian. In mathemati- 
cal terms, this type of inconsistency can occur when the bijective correspondence between 
p-cells of the primal lattice and (n — p)- c ells of the dual lattice (an expression of Poincare 



duality at the level of cellular homology [1051 ] . up to boundary terms) is violated. This is 



typified by 'nonreciprocal' constructions of derivative operators, where the boundary oper- 
ator effecting the spatial derivation on the primal lattice K is not the dual adjoint (or the 
incidence matrix transpose) of the boundary operator on the dual lattice /C: for example, 
the identity Cfj = C^ p (up to boundary terms) used to obtain eq. ( |32l) is violated. One 
basic consequence of this violation is that the resulting discrete equations break time-reversal 
symmetry. Consequently, the numerical solutions will violate energy conservation and pro- 
duce either artificial dissipation or late-time instabilities [22]. Many algorithms developed 
over the years for hyperbolic partial differential equations do indeed viola t e the se proper- 



ties: they are dissipative and cannot be used for long integration times 106| , 107| | . It should 
be noted at this point that lattice field theories invariably break Lorentz covariance and 
many of the continuum Lagrangian symmetries and, as a result, violate conservation laws 
(currents) by virtue of Noether's theorem. For example, angular momentum conservation 
does not hold exactly on the lattice because of the lack of continuous rotational symmetry 
(note that discrete rotational symmetries can still be present). However, this latter type 
of symmetry breaking is of a fundamentally different nature because it is 'controllable', i.e. 
their effect on the computed solutions is made arbitrarily small in the continuum limit. 
More importantly, discrete transcription s of the Noether's theorem can be constructed for 



Lagrangian symmetries on a lattice 13J, 1081 ] . to yield exact conservation laws of (properly 



defined) quantities such as discrete energy and discrete momentum 

(c) Hodge-star inconsistencies: In the third type of inconsistency, we include those that arise 
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in connection with metric proper ties o: 
in the Hodge-star operators |22]. |l09] . 



' the lattice. Because the metric is entirely encoded 



42], such inconsistencies can be simply understood 



as inconsistencies on the construction of discrete Hodge-star operators (or their procedural 
analogues). For example, it is not uncommon for naive disc retiz a tion s in irregular lattices 
to yield asymmetric discrete Hodge operators, as noted in 110] . 111] . Even if symmetry 
is observed, non positive definiteness might en sue that is often associated with portions of 
the lattice with highly skewed or obtuse cells 112] . Lack of either of these properties lead 
to unconditional instabilities that destroy marching-on-time solutions 22]. When very long 
integration times are needed, asymmetry in the discrete Hodge matrices can be a problem 
even if produced at the level of machine rounding-off errors. 



APPENDIX E: Overview of related discretization approaches 

We outline below some discretization programs that rely, one way or another, on tenets 
exposed above. This delineation is mostly informed mostly by applications related to elec- 
trodynamics and not too sharp as the programs share much in common. 

(a) Finite- difference time- domain method: In cubical lattices, the (lowest-order) Whitney 
forms can be represent ed b y means of a product of pulse and 'rooftop' functions on the three 
Cartesian coordinates 1131 ]. This choice, together with the use of low-order quadrature rules 



to compute the Hodge star integrals in ( 1231) and ( |24|) . leads to diagonal matrices [* e ], [* M -i], 
and, consequently, also diagonal [*e]~\ t*^- 1 ] -1 anc ^ s P arse ["£] so that an ultra-local equation 
results for ( 1331) . In this fashion, one obtains a 'matrix- free' algorithm where no linear algebra 
is needed during a marching-on-time solution for t he fields. This prescription recovers Yee's 



finite-difference time-domain scheme 



51], 114J. Conventional FDTD adopts the simplest 



explicit, energy-conserving (symplectic) time-discretization for eqs. (1301) and fl33|) . which can 
be constructed by staggering the electric and magnetic fields in time and replacing time 
derivatives by central differences. Staggering in both space an d tim e is consistent with the 
presence of two staggered hypercubical spacetime lattices 48j, 1151 ]. The staggering in time 
also provides a 0(At 2 ) truncation error. 
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is 



(b) Finite integration technique: The finite integration technique (FIT) [116] . [117 ] 
closely related to FDTD, the main distinction being that, assuming piecewise constant fields 
over each cell, the latter is equivalent to applying the (discrete version) of the generalized 
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Stokes' theorem to the cochains in (|30|) and f|3Tl) . Another difference is that the incidence 
matrices and material (Hodge star) matrices are treated separately in FIT, in a manner 
akin to that exposed in Sections III and IV. Like FDTD, FIT is based on dual staggered 
lattices and, for cubical lattices, it turns out that the lowest-order numerical implementation 
of FIT is equivalent to the lowest-order FDTD. The spatial operators in FIT can all be 
viewed as discrete incarnations of the exterior derivative for the various p, and as such, the 
exact sequence property of the underlying de Rham complex is automatically enforced by 
construction 55|. Historically, F IT g eneralizations to ir regu lar lattices have relied on the 
use of either projection operators 112| or Whitney forms 1191 ] to construct discrete versions 



of the Hodge star operators (or their procedural equivalents); however, these generalizations 
do not necessarily recover the specific form of the discrete Hodge matrix elements expressed 
in (El and (|24j). 

es orig- 



(c) Cell method: Another related discretization pro gram , based on general pr i ncip 



122j.[123|.|l24).|l25|. Even 



inally put forth in |48| , |49[ , |47| , is the Cell method |120j.|l21|. 
though this program does not rely on Whitney forms for constructing discrete Hodge star 
operators (other geometrically-based constructions are used instead), it is nevertheless still 
based upon the use of 'domain-integrated' discrete variables that conform to the notion of 
discrete differential forms or cochains of various degrees and, as such, it is naturally suited for 
irregular lattices. The Cell method also employs metric-free discrete operators that satisfy 
the exactness property of the de Rham complex and make explicit use of a dual lattice (but 
not necessarily barycentric) motivated by the notion of inner and outer orientations. The 
relationships between the various discrete operators and 'domain-integrated' field quantities 
(cochains) in the Cell method are built into general classification diagrams referred to as 
'Tonti diagrams' that reproduce correct commuting diagram properties of the underlying 
operators J47l|.[48|. 

(d) Mimetic finite- differences: 'Mimetic' finite-difference methods, originally developed for 
non-orthogonal hexahedral struct ured l attic e s ('tensor- p roduct gr i ds') and la t er ex t ende d for 
irregular and polyhedral lattices |12jj| , |l27[ , |12s| , |12a| , [l3o| , [13l| , |132j , |133| , 11341 , [13a| also 
share many of the properties exposed above. The thrust here is towards the construction of 
discrete versions of the differential operators divergence, gradient, and curl of vector calcu- 
lus having 'compatible' (in the sense of the exactness property of the underlying de Rham 
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complex) domains and ranges and such that the resulting discrete equations exactly satisfy 
discrete conservation laws. In three dimensions, this naturally leads to the definition of 
three 'natural' operators and three 'adjoint' operators that can be associated with exterior 
derivative d and the codifferential 5, respectively, for p = 1,2,3 (although the exterior cal- 
culus terminology is often not used explicitly in this context). In mimetic finite-differences, 
the discrete analogues of the codifferential operator S are full matrices, and the matrix-free 
character of FDTD is lacking even on orthogonal lattices. 

(e) Compatible discretizations and finite element exterior calculus: In recent years, much 
attention has been devoted to the development of 'compatible discretizations,' an umbrella 
term used to denote spatial discretizations of partial differential equations seeking to provide 
finite element spaces that reproduce the exactness of the u nder l ying de R h am c ompl e x (o r 
the correct cohomology in topologically nontrivial domains) 136| . 1371 ]. 138( |. 1391 ] . 1401 ] . 141 ]. 
In this program, Whitney forms play a role of providing 'conforming' vector-valued func- 
tional (finite element) spaces of Sobolev-type. Specifically, Whitney 1-forms re cove r the 
space of 'Nedelec edge-elements' or curl-conforming Sobolev space H(curl, Q) 1421 ] and 
Whitney 2-forms recover the space of 'Raviart-Thomas elements' or div- conforming Sobolev 
space H(div, fl) 143| . In this regard, a relatively new advance here has been the devel- 
opment of new finite elem ent spaces, beyond those provided by Whitney forms, based on 
the Kos zul complex 1441 ]. The latter is key for the stable discretization of elastodynam- 
ics 1451 ] . Anot her recent approach aimed at the stable discretization of elastodynamics 



is describ ed i n 146] . The link between stability conditions of some mixed finite element 



methods 
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tro dynamics 



an d the complex of Whitney forms 



551] . [147|] . [181] . [191] - J21I] . [231] . [321] . [361] . [611] . [148[| . [149j] . [1501] . [15l|] . 



l as a l ong h istory in t he context of elec 



(f) Discrete exterior calculus: The 'discrete exterior calculus' (DEC) is yet another dis- 
cretization progra m aimed a t developing ah initio consistent discrete models to describe 
field theories jl03| . J152] . J153] . J154J . jl55| . jl56| . This program recognizes the role played by 
discrete differential forms to capture and the need for dual lattices to capture the correct 
physics. Note th at DEC h as focused on the use of a circumcentric dual as opposed to a 
barycentric dual 153], 154] (despite the fact that the former does not admit a metric-free 
construction) and does not emphasize the role of Whitney forms. DEC also recognizes the 
need to address group-valued differential forms, as well as the mathematical objects that 
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exist on the dual-bundle space together with the associated operators (such as contractions 
and Lie derivatives), in connection to discrete problems in mechanics, optimal control, and 



computer vision/graphics [1521 ] . A rec ent d iscussion on obstacles associated with some of 



the DEC underpinnings is provided in 157 ] 
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